library(sf)
Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
library(rgdal)
Please note that rgdal will be retired during October 2023,
plan transition to sf/stars/terra functions using GDAL and PROJ
at your earliest convenience.
See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
rgdal: version: 1.6-7, (SVN revision 1203)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 3.6.2, released 2023/01/02
Path to GDAL shared files: C:/Users/jcolpitt/AppData/Local/R/win-library/4.3/rgdal/gdal
GDAL does not use iconv for recoding strings.
GDAL binary built with GEOS: TRUE
Loaded PROJ runtime: Rel. 9.2.0, March 1st, 2023, [PJ_VERSION: 920]
Path to PROJ shared files: C:\Program Files\PostgreSQL\15\share\contrib\postgis-3.3\proj
PROJ CDN enabled: FALSE
Linking to sp version:1.6-1
To mute warnings of possible GDAL/OSR exportToProj4() degradation,
use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
library(ggplot2)
library(dplyr)
Attaching package: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
library(tidyr)
library(scales)
library(RColorBrewer)
library(units)
udunits database from C:/Users/jcolpitt/AppData/Local/R/win-library/4.3/units/share/udunits/udunits2.xml
library(cowplot)
library(here)
here() starts at C:/Users/jcolpitt/SpatialDataScience/Rprojects/R
okcounty <- st_read(here("data", "ok_counties.shp"), quiet = TRUE)
tpoint <- st_read(here("data", "ok_tornado_point.shp"), quiet = TRUE)
tpath <- st_read(here("data", "ok_tornado_path.shp"), quiet = TRUE)
class(okcounty)
[1] "sf" "data.frame"
glimpse(okcounty)
Rows: 77
Columns: 8
$ STATEFP <chr> "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", …
$ COUNTYFP <chr> "077", "025", "011", "107", "105", "153", "001", "053", "059", "073", "129", "015", "097", "033", "147", "041", "109", "057", "133", "093", "031", "063", "075", "099", "141", "083", "069…
$ COUNTYNS <chr> "01101826", "01101800", "01101793", "01101841", "01101840", "01101864", "01101788", "01101814", "01101817", "01101824", "01101852", "01101795", "01101833", "01101804", "01101861", "01101…
$ AFFGEOID <chr> "0500000US40077", "0500000US40025", "0500000US40011", "0500000US40107", "0500000US40105", "0500000US40153", "0500000US40001", "0500000US40053", "0500000US40059", "0500000US40073", "05000…
$ GEOID <chr> "40077", "40025", "40011", "40107", "40105", "40153", "40001", "40053", "40059", "40073", "40129", "40015", "40097", "40033", "40147", "40041", "40109", "40057", "40133", "40093", "40031…
$ NAME <chr> "Latimer", "Cimarron", "Blaine", "Okfuskee", "Nowata", "Woodward", "Adair", "Grant", "Harper", "Kingfisher", "Roger Mills", "Caddo", "Mayes", "Cotton", "Washington", "Delaware", "Oklahom…
$ LSAD <chr> "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", …
$ geometry <POLYGON [°]> POLYGON ((-95.50766 35.0292..., POLYGON ((-103.0025 36.6751..., POLYGON ((-98.6369 36.16489..., POLYGON ((-96.62486 35.4627..., POLYGON ((-95.80982 36.9419..., POLYGON ((-99.6055…
okcounty_sp <- as(okcounty, 'Spatial')
class(okcounty_sp)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"
okcounty_sf <- st_as_sf(okcounty_sp)
class(okcounty_sf)
[1] "sf" "data.frame"
ggplot(data = okcounty) +
geom_sf(fill = NA)

names(tpoint)
[1] "om" "yr" "mo" "dy" "date" "time" "tz" "st" "stf" "stn"
[11] "mag" "inj" "fat" "loss" "closs" "slat" "slon" "elat" "elon" "len"
[21] "wid" "fc" "geometry"
tpoint_16_21 <- tpoint %>%
filter(yr >= 2016 & yr <= 2021) %>%
select(om, yr, date)
tpath_16_21 <- tpath %>%
filter(yr >= 2016 & yr <= 2021) %>%
select(om, yr, date)
ggplot() +
geom_sf(data = okcounty, fill = NA) +
geom_sf(data = tpoint_16_21) +
theme_bw()

ggplot() +
geom_sf(data = okcounty, fill = NA) +
geom_sf(data = tpath_16_21, color = "red") +
theme_void()

ggplot() +
geom_sf(data = tpoint_16_21,
aes(color = as.factor(yr))) +
geom_sf(data = okcounty, fill = NA) +
scale_color_discrete(name = "Year") +
coord_sf(datum = NA) +
theme_void()

ggplot() +
geom_sf(data = okcounty,
fill = NA,
color = "gray") +
geom_sf(data = tpoint_16_21,size = 0.75) +
facet_wrap(vars(yr), ncol = 2) +
coord_sf(datum = NA) +
theme_void()

countypnt <- st_join(tpoint_16_21, okcounty)
glimpse(countypnt)
Rows: 434
Columns: 11
$ om <dbl> 613662, 613675, 613676, 613677, 613678, 613727, 613728, 613729, 613730, 613731, 613763, 613764, 61376…
$ yr <dbl> 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016,…
$ date <chr> "2016-03-23", "2016-03-30", "2016-03-30", "2016-03-30", "2016-03-30", "2016-04-15", "2016-04-15", "20…
$ STATEFP <chr> "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40",…
$ COUNTYFP <chr> "001", "113", "105", "131", "035", "139", "139", "139", "139", "139", "017", "109", "109", "113", "12…
$ COUNTYNS <chr> "01101788", "01101844", "01101840", "01101853", "01101805", "01101857", "01101857", "01101857", "0110…
$ AFFGEOID <chr> "0500000US40001", "0500000US40113", "0500000US40105", "0500000US40131", "0500000US40035", "0500000US4…
$ GEOID <chr> "40001", "40113", "40105", "40131", "40035", "40139", "40139", "40139", "40139", "40139", "40017", "4…
$ NAME <chr> "Adair", "Osage", "Nowata", "Rogers", "Craig", "Texas", "Texas", "Texas", "Texas", "Texas", "Canadian…
$ LSAD <chr> "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06",…
$ geometry <POINT [°]> POINT (-94.5042 35.6916), POINT (-96.0151 36.2151), POINT (-95.5523 36.7291), POINT (-95.6384 3…
# convert to a non spatial dataframe
countypnt <- st_drop_geometry(countypnt)
countysum <- countypnt %>%
group_by(GEOID) %>%
summarize(tcnt = n())
glimpse(countysum)
Rows: 75
Columns: 2
$ GEOID <chr> "40001", "40005", "40007", "40009", "40011", "40013", "40015", "40017", "40019", "40021", "40023", "4002…
$ tcnt <int> 6, 3, 4, 8, 1, 4, 10, 5, 7, 5, 3, 12, 10, 5, 5, 1, 7, 9, 7, 8, 2, 2, 2, 3, 9, 5, 4, 4, 2, 1, 4, 3, 9, 1,…
countymap <- okcounty %>%
left_join(countysum, by = "GEOID") %>%
replace(is.na(.), 0) %>%
mutate(area = st_area(okcounty),
tdens = 10^6 * 10^3 * tcnt / area) %>%
drop_units()
glimpse(countymap)
Rows: 77
Columns: 11
$ STATEFP <chr> "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40",…
$ COUNTYFP <chr> "077", "025", "011", "107", "105", "153", "001", "053", "059", "073", "129", "015", "097", "033", "14…
$ COUNTYNS <chr> "01101826", "01101800", "01101793", "01101841", "01101840", "01101864", "01101788", "01101814", "0110…
$ AFFGEOID <chr> "0500000US40077", "0500000US40025", "0500000US40011", "0500000US40107", "0500000US40105", "0500000US4…
$ GEOID <chr> "40077", "40025", "40011", "40107", "40105", "40153", "40001", "40053", "40059", "40073", "40129", "4…
$ NAME <chr> "Latimer", "Cimarron", "Blaine", "Okfuskee", "Nowata", "Woodward", "Adair", "Grant", "Harper", "Kingf…
$ LSAD <chr> "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06",…
$ tcnt <dbl> 1, 12, 1, 10, 6, 2, 6, 0, 4, 9, 3, 10, 12, 1, 2, 8, 13, 4, 5, 5, 5, 1, 9, 3, 16, 6, 9, 1, 8, 5, 4, 7,…
$ geometry <POLYGON [°]> POLYGON ((-95.50766 35.0292..., POLYGON ((-103.0025 36.6751..., POLYGON ((-98.6369 36.16489..…
$ area <dbl> 1890663261, 4766283042, 2427121029, 1657249513, 1503893122, 3227039659, 1493619771, 2597526414, 26851…
$ tdens <dbl> 0.5289149, 2.5176851, 0.4120108, 6.0340944, 3.9896452, 0.6197631, 4.0170866, 0.0000000, 1.4896951, 3.…
st_write(countymap, dsn = "oktornado.gpkg",
layer = "countysum",
delete_dsn = TRUE)
st_write(countymap, dsn = "oktornado.geojson",
layeroptions = "RFC7946 = YES")
# proud to say the highest tornado density is my home town, Tulsa County
ggplot(data = countymap) +
geom_sf(aes(fill = tdens)) +
theme_void()

st_geometry_type(okcounty, by_geometry = FALSE)
[1] POLYGON
18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT MULTILINESTRING ... TRIANGLE
okcntrd = st_centroid(countymap)
Warning: st_centroid assumes attributes are constant over geometries
st_geometry_type(okcntrd, by_geometry = FALSE)
[1] POINT
18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT MULTILINESTRING ... TRIANGLE
ggplot() +
geom_sf(data = okcntrd, aes(size = tcnt)) +
geom_sf(data = okcounty, fill = NA) +
theme_void()

ggplot(data = countymap) +
geom_sf(aes(fill = tdens)) +
scale_fill_distiller(name = expression("Tornadoes/1000 km"^2),
palette = "YlOrRd",
breaks = pretty_breaks(),
direction = 1) +
theme_void() +
theme(legend.position = "bottom")

numclas <- 4
qbrks <- seq(0, 1, length.out = numclas + 1)
qbrks
[1] 0.00 0.25 0.50 0.75 1.00
countymap <- countymap %>%
mutate(tdens_c1 = cut(tdens,
breaks = quantile(tdens, breaks = qbrks),
include.lowest = T))
ggplot(data = countymap) +
geom_sf(aes(fill = tdens_c1)) +
scale_fill_brewer(name = expression("Tornadoes/1000 km"^2),
palette = "YlOrRd") +
theme_void() +
theme(legend.position = "bottom")

maxcnt <- max(okcntrd$tcnt)
brkpts <- c(0, 2, 5, 10, maxcnt)
okcntrd <- okcntrd %>%
mutate(tcnt_c1 = cut(tcnt,
breaks = brkpts,
include.lowest = T))
ggplot(data = okcntrd) +
geom_sf(aes(size = tcnt_c1)) +
scale_size_discrete(name="Tornadoes") +
geom_sf(data = okcounty, fill = NA) +
theme_void() +
theme(legend.position = "bottom")
Warning: Using size for a discrete variable is not advised.

ggsave("tornado.png",
width = 6,
height = 4,
dpi = 300)
ggsave("tornado2.png",
width = 15,
height = 10,
units = "cm",
dpi = 100)
choropleth <- ggplot(data = countymap) +
geom_sf(aes(fill = tdens_c1)) +
scale_fill_brewer(name="Density",
palette = "YlOrRd") +
theme_void()
gradsymbol <- ggplot(data = okcntrd) +
geom_sf(aes(size = tcnt_c1)) +
scale_size_discrete(name="Count") +
geom_sf(data = okcounty, fill = NA) +
theme_void()
Warning: Using size for a discrete variable is not advised.
ggsave("choropleth.tiff",
plot = choropleth,
width = 6,
height = 4,
dpi = 300,
compression = "lzw")
ggsave("gradsymbol.tiff",
plot = gradsymbol,
width = 6,
height = 4,
dpi = 300,
compression = "lzw")
plot_grid(choropleth, gradsymbol,
labels = c("A) Choropleth Map",
"B) Graduated Symbol Map",
label_size = 12),
ncol = 1,
hjust = 0,
label_x = 0)

tpath_years <- ggplot() +
geom_sf(data = tpath_16_21,
aes(color = as.factor(yr))) +
geom_sf(data = okcounty, fill = NA) +
scale_color_discrete(name = "Year") +
coord_sf(datum = NA) +
theme_void()
tpoint_years <- ggplot() +
geom_sf(data = tpoint_16_21,
aes(color = as.factor(yr))) +
geom_sf(data = okcounty, fill = NA) +
scale_color_discrete(name = "Year") +
coord_sf(datum = NA) +
theme_void()
ggsave("tpath_years.tiff",
plot = choropleth,
width = 6,
height = 6,
dpi = 300,
compression = "lzw")
ggsave("tpoint_years.tiff",
plot = choropleth,
width = 6,
height = 6,
dpi = 300,
compression = "lzw")
plot_grid(tpath_years, tpoint_years,
labels = c("A) Tornado Paths",
"B) Tornado Points",
label_size = 12),
ncol = 1,
rel_heights = c(1, 1),
align = "hv",
vjust = 5)
ggplot(data = countymap) +
geom_sf(aes(fill = tdens)) +
scale_fill_distiller(name = expression("Tornadoes/1000 km"^2),
palette = "YlOrRd",
breaks = pretty_breaks(),
direction = 1) +
theme_void() +
theme(legend.position = "bottom")

df_join <- countymap %>%
left_join(countypnt, by = "GEOID") %>%
replace(is.na(.), 0) %>%
filter(yr != 0)
ggplot(data = df_join) +
geom_sf(aes(fill = tdens)) +
geom_sf(data = okcounty,
fill = NA,
color = "gray") +
scale_fill_distiller(name = expression("Tornadoes/1000 km"^2),
palette = "YlOrRd",
breaks = pretty_breaks(),
direction = 1) +
facet_wrap(vars(yr), ncol = 2) +
coord_sf(datum = NA) +
theme_void() +
theme(legend.position = "bottom")

url <- "https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_month.geojson"
earthquakes <- readOGR(url)
Warning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among others
OGR data source with driver: GeoJSON
Source: "https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_month.geojson", layer: "2.5_month"
with 1947 features
It has 27 fields
Integer64 fields read as strings: time updated
eqsf <- st_as_sf(earthquakes)
faults_url <- "https://services.arcgis.com/jIL9msH9OI208GCb/ArcGIS/rest/services/Active_Faults/FeatureServer/0/query?where=1%3D1&objectIds=&time=&geometry=&geometryType=esriGeometryEnvelope&inSR=&spatialRel=esriSpatialRelIntersects&resultType=none&distance=0.0&units=esriSRUnit_Meter&relationParam=&returnGeodetic=false&outFields=*&returnGeometry=true&featureEncoding=esriDefault&multipatchOption=xyFootprint&maxAllowableOffset=&geometryPrecision=&outSR=&defaultSR=&datumTransformation=&applyVCSProjection=false&returnIdsOnly=false&returnUniqueIdsOnly=false&returnCountOnly=false&returnExtentOnly=false&returnQueryGeometry=false&returnDistinctValues=false&cacheHint=false&orderByFields=&groupByFieldsForStatistics=&outStatistics=&having=&resultOffset=&resultRecordCount=&returnZ=false&returnM=false&returnExceededLimitFeatures=true&quantizationParameters=&sqlFormat=none&f=pgeojson&token="
faults <-readOGR(faults_url)
Warning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among others
OGR data source with driver: GeoJSON
Source: "https://services.arcgis.com/jIL9msH9OI208GCb/ArcGIS/rest/services/Active_Faults/FeatureServer/0/query?where=1%3D1&objectIds=&time=&geometry=&geometryType=esriGeometryEnvelope&inSR=&spatialRel=esriSpatialRelIntersects&resultType=none&distance=0.0&units=esriSRUnit_Meter&relationParam=&returnGeodetic=false&outFields=*&returnGeometry=true&featureEncoding=esriDefault&multipatchOption=xyFootprint&maxAllowableOffset=&geometryPrecision=&outSR=&defaultSR=&datumTransformation=&applyVCSProjection=false&returnIdsOnly=false&returnUniqueIdsOnly=false&returnCountOnly=false&returnExtentOnly=false&returnQueryGeometry=false&returnDistinctValues=false&cacheHint=false&orderByFields=&groupByFieldsForStatistics=&outStatistics=&having=&resultOffset=&resultRecordCount=&returnZ=false&returnM=false&returnExceededLimitFeatures=true&quantizationParameters=&sqlFormat=none&f=pgeojson&token=", layer: "OGRGeoJSON"
with 16121 features
It has 6 fields
faults_sf <- st_as_sf(faults)
pal <- colorBin(
palette = "Spectral",
domain = eqsf$mag,
reverse = TRUE,
bins = 5
)
leaflet() %>%
addTiles() %>%
setView(-117.841293, 46.195042, 3) %>%
addProviderTiles(providers$CartoDB, group = "Grayscale") %>%
addProviderTiles(providers$Esri.WorldTerrain, group = "Terrain") %>%
addPolylines(data = faults_sf, group = "vectorData") %>%
addCircleMarkers(data = eqsf,
fillColor = ~pal(mag),
radius = ~eqsf$mag * 2,
stroke = FALSE,
color = "Spectral",
fillOpacity = 0.6,
popup = paste0(
"<strong>Title:</strong> ", eqsf$title,
"<br><strong>Magnitude:</strong> ", eqsf$mag,
"<br><strong>Intensity:</strong> ", eqsf$mmi,
"<br><strong>Sig:</strong> ", eqsf$sig
),
group = "vectorData"
) %>%
addLegend(pal = pal, values = eqsf$mag, position = "bottomleft", title = "Magnitude") %>%
addLayersControl(overlayGroups = c("vectorData"), baseGroups = c("Terrain", "Grayscale"))